Comparaison du MIC de matrices de spectres Données human Serum prétraitées par SOAP par deux méthodes:

  1. Full : Peps avec l’ensemble des étapes
  2. Manuel : étapes appliquées manuellement avec Topspin par l’équipe ULG

1 Lecture des matrices des spectres et initiatilisation des variables utiles

#  Choix du répertoire où sont les matrices de données
pathdir="~/Dropbox/Partage_OMICS_COURSE/Data_stat2340/Human_Serum"
#  Definition des matrices à lire dans ce répertoire
data_used=c("HumSerum_full_PEPS.csv","HumSerum_Manual_500P.csv")
spectral_matrices_names=c("HumSerum_full_PEPS","HumSerum_Manual_500P")
nmatrices=length(data_used)
# Lecture des matrices de données.  
## On les stocke dans spectral_matrices_list
## On remplace les noms des colonnes par des noms sans le X
spectral_matrices_list=list() 
for(i in 1:nmatrices)
{
spectral_matrices_list[[spectral_matrices_names[i]]]=as.matrix(read.table(file.path(pathdir,data_used[i]),header=T,sep=";",row.names=1,dec="."))
dimnames(spectral_matrices_list[[i]])[[2]]=as.numeric(substring(dimnames(spectral_matrices_list[[i]])[[2]],first=2))
cat("\n Lecture de la matrice : ",spectral_matrices_names[i])
}
## 
##  Lecture de la matrice :  HumSerum_full_PEPS
##  Lecture de la matrice :  HumSerum_Manual_500P
# On prend la matrice des spectres communs dans les 2 bases (si ce n'est pas le cas au départ), on eleve les blancs des noms 
dimnames(spectral_matrices_list[[1]])[[1]]=gsub(" ","",dimnames(spectral_matrices_list[[1]])[[1]]) # on enleve des blancs des noms des spectres.  
namespectra=dimnames(spectral_matrices_list[[1]])[[1]]
countspectra=matrix(nrow=nmatrices+1,ncol=2)
dimnames(countspectra)=list(c(spectral_matrices_names,"Common Spectra"),c("Nb of spectra","Nb Buckets"))
countspectra[1,1]=length(namespectra)
countspectra[1,2]=dim(spectral_matrices_list[[1]])[2]
for(i in 2:nmatrices)
{
dimnames(spectral_matrices_list[[i]])[[1]]=gsub(" ","",dimnames(spectral_matrices_list[[i]])[[1]]) 
newnames=dimnames(spectral_matrices_list[[i]])[[1]]
compnames=is.element(namespectra,newnames) 
countspectra[i,1]=length(newnames)
countspectra[i,2]=dim(spectral_matrices_list[[i]])[2]
namespectra=namespectra[compnames]
}  
countspectra[nmatrices+1,1]=length(namespectra)
pander(countspectra)
  Nb of spectra Nb Buckets
HumSerum_full_PEPS 32 500
HumSerum_Manual_500P 31 500
Common Spectra 31 NA
# On extrait les bons spectres et on les met dans le même ordre
for(i in 1:nmatrices)
{
spectral_matrices_list[[i]]=spectral_matrices_list[[i]][namespectra,]
} 
# Creation des variables qui définissent les groupes
# Pour les données human serum c'est le 5ième charactère du nom du spectre
groupes=as.numeric(substring(dimnames(spectral_matrices_list[[1]])[[1]], 5, 5))
names_spectra=substring(dimnames(spectral_matrices_list[[1]])[[1]],1,5)
repetition=as.numeric(substring(dimnames(spectral_matrices_list[[1]])[[1]], 2, 2))
# 3. on dessine les spectres
for(i in 1:nmatrices)
{
print(spectral_matrices_names[i])  
Draw(spectral_matrices_list[[i]],type.draw="signal",subtype="stacked",num.stacked=4)
}
## [1] "HumSerum_full_PEPS"

## [1] "HumSerum_Manual_500P"

2 Etude de la répétabilité spectrale

2.1 Between and Within group inertia

# Impression des résultats des calculs d'inertie
cat("Inerties non standardisées")
## Inerties non standardisées
pander(list(Inerties_non_standardisees=Inertiavalues,Inerties_en_PC=Inertiapc,Inerties_within_groups=InertiaPerGroup))
  • Inerties_non_standardisees:

      BI WI TI
    HumSerum_full_PEPS 6366 702.3 7069
    HumSerum_Manual_500P 0.02525 0.09521 0.1205
  • Inerties_en_PC:

      BI WI TI
    HumSerum_full_PEPS 90.06 9.936 100
    HumSerum_Manual_500P 20.96 79.04 100
  • Inerties_within_groups:

      HumSerum_full_PEPS HumSerum_Manual_500P
    1 168.8 0.02751
    2 256.4 0.0251
    3 84.12 0.02492
    4 192.9 0.01768
    Total 702.3 0.09521

2.2 PCA

2.2.1 Eigenvalues

2.2.2 PCA Score Plots

draw.scores12 = vector("list")
draw.scores34 = vector("list")
for(i in 1:nmatrices)
{
draw.scores12[[spectral_matrices_names[i]]]=MBXUCL::DrawScores(PCA.res[[i]], drawNames=TRUE, type.obj = "PCA",createWindow=FALSE, main = paste0("PCA score plot for ",spectral_matrices_names[i]),pch = groupes, col = groupes,axes =c(1,2))
draw.scores34[[spectral_matrices_names[i]]]=MBXUCL::DrawScores(PCA.res[[i]], drawNames=TRUE, type.obj = "PCA",createWindow=FALSE, main = paste0("PCA score plot for ",spectral_matrices_names[i]),pch = groupes,col = groupes, axes =c(3,4))
}
draw.scores12
## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

draw.scores34
## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

2.2.3 PCA Loading Plots

draw.loadings12 = vector("list")
draw.loadings34 = vector("list")
for(i in 1:nmatrices)
{
draw.loadings12[[spectral_matrices_names[i]]] = MBXUCL::DrawLoadings(PCA.res[[i]],  type.obj = "PCA",createWindow=FALSE, main = paste0("PCA loadings plot for", spectral_matrices_names[i]),axes = c(1:2), loadingstype="l", num.stacked = 2)[[1]]
draw.loadings34[[spectral_matrices_names[i]]] = MBXUCL::DrawLoadings(PCA.res[[i]],  type.obj = "PCA",createWindow=FALSE, main = paste0("PCA loadings plot for", spectral_matrices_names[i]),axes = c(3:4), loadingstype="l", num.stacked = 2)[[1]]
}
draw.loadings12
## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

draw.loadings34
## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

2.3 Unsupervised clustering (MIC)

nClust = length(unique(groupes)) # nombre de clusters à rechercher
clustres=matrix(0,nrow=nmatrices,ncol=8)
dimnames(clustres)=list(spectral_matrices_names,c("DunnW","DunnKM","DBW","DBKM","RandW","RandKM","AdjRandW","AdjRandKM"))
for(i in 1:nmatrices)
{
print(spectral_matrices_names[i]) 
ClustMIC.res = MBXUCL::ClustMIC(Intensities = spectral_matrices_list[[i]], nClust = nClust, Trcl = groupes, Dendr = TRUE)
clustres[i,]=as.numeric(ClustMIC.res[1,1:8])
}
## [1] "HumSerum_full_PEPS"

## [1] "HumSerum_Manual_500P"

pander(clustres)
Table continues below
  DunnW DunnKM DBW DBKM RandW RandKM
HumSerum_full_PEPS 0.8246 0.8246 0.6776 0.6776 1 1
HumSerum_Manual_500P 0.9332 0.9332 0.7824 0.7824 0.7054 0.7054
  AdjRandW AdjRandKM
HumSerum_full_PEPS 1 1
HumSerum_Manual_500P 0.2324 0.2324

2.4 PLS-DA

nLVPLSDA = 4  # nombre du composantes du PLSDA
nrep=nlevels(as.factor(groupes)) # Nombre de réponses
PLSDA.res=list()
perf.plsda.RMSEP=matrix(nrow=nmatrices,ncol=nrep)
dimnames(perf.plsda.RMSEP)=list(spectral_matrices_names,paste0("Y",1:nLVPLSDA))
perf.plsda.R2=perf.plsda.RMSEP
perf.plsda.Q2=perf.plsda.RMSEP
for(i in 1:nmatrices)
{
cat("PLS-DA pour ",spectral_matrices_names[i])  
PLSDA.res [[i]]= PLSDA(x = spectral_matrices_list[[i]], y = groupes, nLV = nLVPLSDA, drawRMSEP = TRUE)
perf.plsda.RMSEP[i,] = PLSDA.res[[i]]$RMSEP
perf.plsda.R2[i,]=PLSDA.res[[i]]$R2
perf.plsda.Q2[i,]=PLSDA.res[[i]]$Q2
}
## PLS-DA pour  HumSerum_full_PEPS

## PLS-DA pour  HumSerum_Manual_500P

pander(list(RMSEP=perf.plsda.RMSEP,R2=perf.plsda.R2,Q2=perf.plsda.Q2))
  • RMSEP:

      Y1 Y2 Y3 Y4
    HumSerum_full_PEPS 0.1145 0.0792 0.1061 0.1117
    HumSerum_Manual_500P 0.3953 0.2478 0.426 0.2985
  • R2:

      Y1 Y2 Y3 Y4
    HumSerum_full_PEPS 0.9728 0.9873 0.971 0.9662
    HumSerum_Manual_500P 0.619 0.8414 0.6475 0.6968
  • Q2:

      Y1 Y2 Y3 Y4
    HumSerum_full_PEPS 0.9316 0.9672 0.9356 0.9349
    HumSerum_Manual_500P 0.1839 0.6794 -0.0381 0.5347

2.4.1 Score plots de la PLS_DA

PLSDA.scores= vector("list")
for(i in 1:nmatrices)
{ 
PLSDA.scores[[spectral_matrices_names[i]]]=DrawScores(PLSDA.res[[i]], drawNames=TRUE, type.obj = c("PLSDA"),
        createWindow=FALSE, main = paste0("PLSDA score plot for ", spectral_matrices_names[i]),pch = groupes,col = groupes, axes =c(1,2))
}
PLSDA.scores
## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

2.4.2 Loading plots de la PLS_DA

PLSDA.loadings= vector("list")
for(i in 1:nmatrices)
{ 
PLSDA.loadings[[spectral_matrices_names[i]]]=DrawLoadings(PLSDA.res[[i]],  type.obj = "PLSDA",
        createWindow=FALSE, main = paste0("PLSDA loadings plot for", spectral_matrices_names[i]),
        axes = c(1:2),  loadingstype="l", num.stacked = 2)
}
PLSDA.loadings
## $HumSerum_full_PEPS
## $HumSerum_full_PEPS[[1]]

## 
## 
## $HumSerum_Manual_500P
## $HumSerum_Manual_500P[[1]]